######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate Figure B.1 - Figure B.3
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#Load bpCausal package
library(bpCausal)

#Other packages used
library(tidyverse)
library(RColorBrewer)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

###################################
#Estimation for Figure B.1 - B.3
###################################
#----------------
#(1) Load data
#----------------
load("data/data_main.RData")

#----------------
#(2) Setup
#----------------
#Setup
df.use = df.main
covs.all = c("gdppc_log_lag",
             "population_log_lag",
             "election_either_lag",
             "fdi_gdp_lag",
             "debt_gni_lag",
             "oda_gni_lag",
             "polity2_lag",
             "unsc_lag",
             "IdealPointDistance_lag")
X.use = Z.use = A.use = covs.all
dv.list = c("hard_project_count_all")

#----------------
#(3) Estimation: Placebo 2014
#~3 minutes
#----------------
#Create AIIB Founder * 2014
df.use$aiib_founder_2014 = ifelse(df.use$aiib_founder == 1 & df.use$year >= 2014, 1, 0)
D.use = "aiib_founder_2014"

result.b1 = bpcausal.group(dv.list = dv.list,
                           df.use = df.use,
                           D.use = D.use,
                           X.use = X.use,
                           Z.use = Z.use,
                           A.use = A.use)

#----------------
#(4) Estimation: Placebo 2012
#~3 minutes
#----------------
#Create AIIB Founder * 2012
df.use$aiib_founder_2012 = ifelse(df.use$aiib_founder == 1 & df.use$year >= 2012, 1, 0)
D.use = "aiib_founder_2012"

result.b2 = bpcausal.group(dv.list = dv.list,
                           df.use = df.use,
                           D.use = D.use,
                           X.use = X.use,
                           Z.use = Z.use,
                           A.use = A.use)

#----------------
#(5) Estimation: Placebo 2010
#~3 minutes
#----------------
#Create AIIB Founder * 2010
df.use$aiib_founder_2010 = ifelse(df.use$aiib_founder == 1 & df.use$year >= 2010, 1, 0)
D.use = "aiib_founder_2010"

result.b3 = bpcausal.group(dv.list = dv.list,
                           df.use = df.use,
                           D.use = D.use,
                           X.use = X.use,
                           Z.use = Z.use,
                           A.use = A.use)

###################################
#Produce Figure B.1
###################################
#-------------
#(1) Prepare results
#-------------
result.14.use = effSummary(result.b1$hard_project_count_all)

#-------------
#(2) Parameters
#-------------
#Titles
xtitle = ""
ytext = "Effect on New World Bank Infrastructure Projects"

#Limits
ylim = c(-3, 3)
newx = 1992:2019
t.time = 2016

#Cex
cex = 2

#-------------
#(3) Plot
#-------------
#Setup
png(filename = "figure/Figure_B1.png",
    height = 800,
    width = 1200)
par(mar = c(2, 5, 5, 1),
    lend = 1)

#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = "dotted",
       lwd = cex)

abline(v = 2014,
       col = "gray",
       lty = "dashed",
       lwd = cex)

abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = result.14.use$est.eff$estimated_ATT,
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(result.14.use$est.eff$estimated_ATT_ci_l), 
              result.14.use$est.eff$estimated_ATT_ci_u),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Credible Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")

#Close
dev.off()

###################################
#Produce Figure B.2
###################################
#-------------
#(1) Prepare results
#-------------
result.12.use = effSummary(result.b2$hard_project_count_all)

#-------------
#(2) Parameters
#-------------
#Titles
xtitle = ""
ytext = "Effect on New World Bank Infrastructure Projects"

#Limits
ylim = c(-3, 3)
newx = 1992:2019
t.time = 2016

#Cex
cex = 2

#-------------
#(3) Plot
#-------------
#Setup
png(filename = "figure/Figure_B2.png",
    height = 800,
    width = 1200)
par(mar = c(2, 5, 5, 1),
    lend = 1)

#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = "dotted",
       lwd = cex)

abline(v = 2012,
       col = "gray",
       lty = "dashed",
       lwd = cex)

abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = result.12.use$est.eff$estimated_ATT,
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(result.12.use$est.eff$estimated_ATT_ci_l), 
              result.12.use$est.eff$estimated_ATT_ci_u),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Credible Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")

#Close
dev.off()

###################################
#Produce Figure B.3
###################################
#-------------
#(1) Prepare results
#-------------
result.10.use = effSummary(result.b3$hard_project_count_all)

#-------------
#(2) Parameters
#-------------
#Titles
xtitle = ""
ytext = "Effect on New World Bank Infrastructure Projects"

#Limits
ylim = c(-3, 3)
newx = 1992:2019
t.time = 2016

#Cex
cex = 2

#-------------
#(3) Plot
#-------------
#Setup
png(filename = "figure/Figure_B3.png",
    height = 800,
    width = 1200)
par(mar = c(2, 5, 5, 1),
    lend = 1)

#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = "dotted",
       lwd = cex)

abline(v = 2010,
       col = "gray",
       lty = "dashed",
       lwd = cex)

abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = result.10.use$est.eff$estimated_ATT,
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(result.10.use$est.eff$estimated_ATT_ci_l), 
              result.10.use$est.eff$estimated_ATT_ci_u),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Credible Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")

#Close
dev.off()
